---
title: "Supplementary Materials for A Comment on 'Instrumentally Inclusive: The Political Psychology of Homonationalism' (Turnbull-Dugarte and López Ortega, 2024)."
author: 
  - name: "Daniel de Kadt"
    affiliation: "London School of Economics and Political Science"
    email: "d.n.de-kadt@lse.ac.uk"
citeproc: true
bibliography: references.bib
crossref:
  custom:
    - kind: float
      key: smfig
      latex-env: smfig
      reference-prefix: Figure SM
      space-before-numbering: false
      latex-list-of-description: Appendix Figure
    - kind: float
      key: smtbl
      latex-env: smtbl
      reference-prefix: Table SM 
      space-before-numbering: false
      latex-list-of-description: Appendix Table
format: 
  typst:
    toc: true
    mainfont: "Computer Modern"


---

```{r setup, echo = FALSE, include = FALSE, warnings = FALSE}
# Note to replicators from DdK: As much as possible, I have attempted throughout to avoid adjusting the original authors' code. Most changes are documented as comments. Please run this chunk once to install any missing packages using groundhog. 
# load libraries
library(groundhog)

pkgs <- c("tidyverse", "jtools", "ggpubr", "ggrepel", "patchwork", "gt", "modelsummary", "interactions", "margins", "skimr", "survey", "estimatr")
groundhog::groundhog.library(pkgs, "2024-06-01")

# modelsummary options
options(modelsummary_model_labels = "Computer Modern")

# set seed exactly as per replication materials
set.seed(1)
```

```{r data ingestion, echo = FALSE, include = FALSE, warnings = FALSE}
# load both datasets
uk <- read_csv("study1_data.csv")
load("study2_data.Rda")

# set color palette as per replication materials
colors<- c("#205C8A", "#d11141")

# cleaning per replication materials
uk <- uk%>% 
  mutate(treat= as.factor(treatment),
         treatnum= as.numeric(treatment),
         gender= as.factor(gender),
         degree= as.factor(degree),
         nonwhite= as.factor(nonwhite),
         queer= as.factor(queer),
         relig= as.factor(relig),
         religion= as.factor(religion),
         race= as.factor(race),
         fourarm= as.factor(fourarm),
         immbelow= as.factor(immbelow),
         imm3= as.factor(imm3),
         region= as.factor(region),
         voterecall= as.factor(voterecall),
         brexit= as.factor(brexit),
         ideology= as.factor(ideology),
         agecat= as.factor(agecat))

# cleaning per replication materials
spain <- spain%>% 
  mutate(treat= as.factor(treat),
         treatnum= as.numeric(treat),
         gender= as.factor(gender),
         supportcat= as.factor(support),
         agecat= as.factor(agecat),
         child= as.factor(child),
         immdum= as.factor(immdum),
         imm5= as.factor(imm5),
         imm3= as.factor(imm3),
         foreignborn= as.factor(foreignborn),
         CCAA= as.factor(CCAA),
         queer= as.factor(queer))
```

```{r uk analysis and subsetting, echo = FALSE, include = FALSE, warnings = FALSE}
# subset uk data per replication materials
treat <- subset(uk, treatnum==1)
control <- subset(uk, treatnum==0)
proimm <- subset(uk, immbelow==0)
noproimm <- subset(uk, immbelow==1)

# analyse uk data per replication materials
modelsub1<- lm(support ~ treat, data=proimm)
proimm$predictedb <- predict(modelsub1, proimm)

modelsub2<- lm(support ~ treat, data=noproimm)
noproimm$predictedb <- predict(modelsub2, noproimm)

# subset again per replication materials
treatsub1 <- subset(proimm, treatnum==1)
controlsub1 <- subset(proimm, treatnum==0)
treatsub2 <- subset(noproimm, treatnum==1)
controlsub2 <- subset(noproimm, treatnum==0)
```

```{r spain analysis and subsetting, echo = FALSE, include = FALSE, warnings = FALSE}
# subset spain data per replication materials
treatES <- subset(spain, treat==1)
controlES <- subset(spain, treat==0)
proimmES <- subset(spain, immdum==1)
noproimmES <- subset(spain, immdum==0)

# analyse spain data, but use lm() not glm() so that Spain and UK analyses are exactly the same:
modelsub1ES <- lm(support ~ treat, weight=nationalweight, data=proimmES)
proimmES$predictedb <- predict(modelsub1ES, proimmES)

modelsub2ES <- lm(support ~ treat, weight=nationalweight, data=noproimmES)
noproimmES$predictedb <- predict(modelsub2ES, noproimmES)

# subset again per replication materials
treatsub1ES <- subset(proimmES, treat==1)
controlsub1ES <- subset(proimmES, treat==0)
treatsub2ES <- subset(noproimmES, treat==1)
controlsub2ES <- subset(noproimmES, treat==0)
```

```{r spain replication no-weights, echo = FALSE, include = FALSE, warnings = FALSE}
# computational reproduction of Table A9:
models_rep <- list(
  'Base' = lm(support ~ treat + imm_1, data=spain, weight = nationalweight),
  'Interaction model' = lm(support ~ treat*imm_1, data=spain, weight = nationalweight),
  'Pro-immigration only' = lm(support ~ treat, data=proimmES, weight = nationalweight),
  'Anti-immigration only' = lm(support ~ treat, data=noproimmES, weight = nationalweight)
)

# reproduction by analysis, with varying researcher choices:
models_base <- list(
  'Replication' = lm(support ~ treat + imm_1, data=spain, weight = nationalweight),
  'No Weights' = lm(support ~ treat + imm_1, data=spain),
  'Weighted Robust SE' = estimatr::lm_robust(support ~ treat + imm_1, data=spain, weight = nationalweight, se_type = "HC3"),
  'Unweighted Robust SE' = estimatr::lm_robust(support ~ treat + imm_1, data=spain, se_type = "HC3"),
  'Weighted Survey-Robust SE' = survey::svyglm(support ~ treat + imm_1, design=svydesign(ids=~1, weights=~nationalweight, data=spain))
)

models_int <- list(
  'Replication'  = lm(support ~ treat*imm_1, data=spain, weight = nationalweight),
  'No Weights' = lm(support ~ treat*imm_1, data=spain),
  'Weighted Robust SE' = estimatr::lm_robust(support ~ treat*imm_1, data=spain, weight = nationalweight, se_type = "HC3"),
  'Unweighted Robust SE' = estimatr::lm_robust(support ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Weighted Survey-Robust SE' = survey::svyglm(support ~ treat*imm_1, design=svydesign(ids=~1, weights=~nationalweight, data=spain))
)

models_proimmes <- list(
  'Replication'  = lm(support ~ treat, data=proimmES, weight = nationalweight),
  'No Weights' =  lm(support ~ treat, data=proimmES),
  'Weighted Robust SE' = estimatr::lm_robust(support ~ treat, data=proimmES, weight = nationalweight, se_type = "HC3"),
  'Unweighted Robust SE' = estimatr::lm_robust(support ~ treat, data=proimmES, se_type = "HC3"),
  'Weighted Survey-Robust SE' = survey::svyglm(support ~ treat, design=svydesign(ids=~1, weights=~nationalweight, data=proimmES))
)

models_noproimmes <- list(
  'Replication'  = lm(support ~ treat, data=noproimmES, weight = nationalweight),
  'No Weights' =  lm(support ~ treat, data=noproimmES),
  'Weighted Robust SE' = estimatr::lm_robust(support ~ treat, data=noproimmES, weight = nationalweight, se_type = "HC3"),
  'Unweighted Robust SE' = estimatr::lm_robust(support ~ treat, data=noproimmES, se_type = "HC3"),
  'Weighted Survey-Robust SE' = survey::svyglm(support ~ treat, design=svydesign(ids=~1, weights=~nationalweight, data=noproimmES))
)

# ancillary mechanism test (A11):
mech <- list(
  'EU norms' = estimatr::lm_robust(pride_valoresUE ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Western liberal values' = estimatr::lm_robust(pride_libertadOCC ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Green politics' = estimatr::lm_robust(pride_verde ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Domestic violence protections' = estimatr::lm_robust(pride_viomach ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Spanish flag' = estimatr::lm_robust(pride_bandera ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Spanish military efforts' = estimatr::lm_robust(pride_mili ~ treat*imm_1, data=spain, se_type = "HC3")
)

```


```{r weight bin prep, echo = FALSE, include = FALSE, warnings = FALSE}
# create bins of weights:
spain <- spain %>%
  mutate(weightbins = case_when(
    nationalweight >= 0 & nationalweight < 0.01 ~ "low",
    nationalweight >= 0.1 & nationalweight < 3 ~ "medium",
    nationalweight >= 3 ~ "high",
    TRUE ~ "other"
  ))

# re-subset
proimmES <- subset(spain, immdum==1)
noproimmES <- subset(spain, immdum==0)

# build tables - anti-immigration subsample
models_hetfx_anti <- list(
  'Low Weights' = estimatr::lm_robust(support ~ treat, data=noproimmES[noproimmES$weightbins=="low",], se_type = "HC3"),
  'Mid Weights' = estimatr::lm_robust(support ~ treat, data=noproimmES[noproimmES$weightbins=="medium",], se_type = "HC3"),
  'High Weights' = estimatr::lm_robust(support ~ treat, data=noproimmES[noproimmES$weightbins=="high",], se_type = "HC3")
)

# build tables - pro-immigration subsample
models_hetfx_pro <- list(
  'Low Weights' = estimatr::lm_robust(support ~ treat, data=proimmES[proimmES$weightbins=="low",], se_type = "HC3"),
  'Mid Weights' = estimatr::lm_robust(support ~ treat, data=proimmES[proimmES$weightbins=="medium",], se_type = "HC3"),
  'High Weights' = estimatr::lm_robust(support ~ treat, data=proimmES[proimmES$weightbins=="high",], se_type = "HC3")
)


# build tables - interaction hetfx
models_hetfx_int <- list(
  'Low Weights' = estimatr::lm_robust(support ~ treat*imm_1, data=spain[spain$weightbins=="low",], se_type = "HC3"),
  'Mid Weights' = estimatr::lm_robust(support ~ treat*imm_1, data=spain[spain$weightbins=="medium",], se_type = "HC3"),
  'High Weights' = estimatr::lm_robust(support ~ treat*imm_1, data=spain[spain$weightbins=="high",], se_type = "HC3")
)
```


```{r uk interaction analysis, echo = FALSE, include = FALSE, warnings = FALSE}
# correct the code to use a linear regression lm() and not logistic regression. 
model1 <- lm(support ~ treat*imm_1, data=uk)
# use lm_robust for robust SEs that can be used by margins(). Set se_type = "HC3" to be consistent with summ(.,robust=TRUE). Point estimates are of course numerically identical bar rounding.
model1_robust <- estimatr::lm_robust(support ~ treat*imm_1, data=uk, se_type = "HC3")

```


```{r spain interaction analysis, echo = FALSE, include = FALSE, warnings = FALSE}
# correct the code to use a linear regression lm() and not logistic regression. 
modelES <- lm(support ~ treat*imm_1, data=spain, weight=nationalweight)
# use lm_robust for robust SEs that can be used by margins(). Set se_type = "HC3" to be consistent with summ(.,robust=TRUE). Point estimates are of course numerically identical bar rounding.
modelES_robust <- estimatr::lm_robust(support ~ treat*imm_1, data=spain, weight=nationalweight, se_type = "HC3")

```

{{< pagebreak >}}


# Note on Computational Reproducibility

This appendix reports on computational reproducibility. The paper is almost perfectly computationally reproducible from the processed data provided. That is, the authors provide the code and data required to (almost exactly) reproduce all results in the paper and supplementary materials, and the code runs (largely) without error. However, the data have evidently been pre-processed (e.g. the creation of the weights which were likely not provided by Prolific), and neither the raw data nor this pre-processing code are provided. I note a few minor issues encountered below:

- The `modelsummary::modelsummary` function does not support the use of `robust = TRUE` as an argument for returning robust standard errors. Despite this, the standard errors in the paper are robust where this argument is used (perhaps a prior version supported this). The correct result should be achieved with `vcov = "HC3"` instead. 

- The `ritest` package not available on CRAN, and must be downloaded using: `remotes::install_github("grantmcdermott/ritest")`. Ideally this would be noted in the materials. 

- The `starbility` package not available on CRAN, and must be downloaded using:
`remotes::install_github('https://github.com/AakaashRao/starbility')`.  Ideally this would be noted in the materials. 

- When using `ggsave()` it would be wise to hard-code the dimensions of the output figure. This function defaults to the current dimensions of the plot window in RStudio (or whatever graphics device is currently active), which will vary by user. 

- In the script `study1_summarystats.R` the code erroneously asks for `UKdata_analysis.csv` which is not provided in the replication archive. The correct file is `study1_data.csv` -- when correcting this the results are reproduced. There are also some slight inconsistencies in the tables produced by this code and the tables in the supplementary materials (e.g. there are rows returns by the code for Table A.3. that are not in the supplementary materials). There is a similar discrepancy in the code that produces Table A.5. in `study2_summarstats.R`. In this file there are two lines that seem to produce summary statistics tables (one that is weighted, one that is not), but only one is in the supplementary materials (the weighted one). This is not indicated in the table.

- For the power analyses, while I was able to reproduce the results in the supplementary information by running the provided code, I was not able to trace the input values myself. These values are hard-coded in the replication code, but it is unclear where the values come from (they do not appear to come directly from the replication data, based on my cursory explorations).

- It is worth noting that the replication materials do not include any data processing code. Though this is not generally required by journals at present and is very rare in the discipline, it would be particularly useful for understanding some of the issues with regards to the weights in the Spain experiment.

{{< pagebreak >}}

# Regression Results and Additional Analyses

This SI includes a series of reproductions and analyses of the authors' data, with a primary focus on study 2. 

## Re-analyses of main specifications

In this subsection I report the primary re-analyses of study 2, varying both the weighting choices and the standard error estimation. In @smtbl-some-table I first reproduce the results from Table A9 in the published paper's supplementary materials. This table reports results for four different specifications -- the base model, the interaction model, the pro-immigration only sample, and the anti-immigration only sample. Results from Table A9 are reproduced exactly.


::: {#smtbl-some-table}

```{r spain replication weight table, echo = FALSE, warnings = FALSE}
#| label: smtbl-some-table

modelsummary(models_rep, output = "gt", 
             coef_map = c('treat1' = 'Treatment', 'imm_1' = 'Immigration', 'treat1:imm_1' = 'Treatment x Immigration', '(Intercept)' = 'Intercept'),
             gof_omit = "BIC|AIC|R2 Adj.|F|RMSE|Log.Lik.",  
             stars = c('*'=.1, "**"=.05, "***"=.01)) 
```

Spain Experiment: Regression Results With Weights and Classical (Non-Robust) Standard Errors

:::

In @smtbl-base-models I vary the weighting and standard error choices for the base model. I do the same for the interaction models in @smtbl-int-models, for the pro-immigration sample in @smtbl-proimmes-models, and finally for the anti-immigration sample in @smtbl-noproimmes-models. In each case I also reproduce the results from Table A9 as the first column ('replication').

::: {#smtbl-base-models}

```{r spain replication base models, echo = FALSE, warnings = FALSE}
#| label: smtbl-base-models

suppressWarnings(
  modelsummary(models_base, output = "gt", 
             coef_map = c('treat1' = 'Treatment', 'imm_1' = 'Immigration', '(Intercept)' = 'Intercept'),
             gof_omit = "BIC|AIC|R2 Adj.|F|RMSE|Log.Lik.",  
             stars = c('*'=.1, "**"=.05, "***"=.01), 
             add_rows = tribble(~term,  ~Base,  ~Base, ~Base, ~Base, ~Base,
               "Weighted", "Yes", "No", "Yes", "No", "Yes", 
               "HC3 Robust SEs", "No", "No", "Yes", "Yes", "No",
               "Survey Robust SEs", "No", "No", "No", "No", "Yes"
             ))
)

```

Spain Experiment: Base Model (Table A9, Column 1) Sensitivity

:::

::: {#smtbl-int-models}

```{r spain replication interaction models, echo = FALSE, warnings = FALSE}
#| label: smtbl-int-models

suppressWarnings(
  modelsummary(models_int, output = "gt", 
             coef_map = c('treat1' = 'Treatment', 'imm_1' = 'Immigration', 'treat1:imm_1' = 'Treatment x Immigration', 
                          '(Intercept)' = 'Intercept'),
             gof_omit = "BIC|AIC|R2 Adj.|F|RMSE|Log.Lik.",  
             stars = c('*'=.1, "**"=.05, "***"=.01), 
             add_rows = tribble(~term,  ~Base,  ~Base, ~Base, ~Base, ~Base,
               "Weighted", "Yes", "No", "Yes", "No", "Yes", 
               "HC3 Robust SEs", "No", "No", "Yes", "Yes", "No",
               "Survey Robust SEs", "No", "No", "No", "No", "Yes"
             ))
)

```

Study 2 (Spain): Interaction Model (Table A9, Column 2) Sensitivity

:::

::: {#smtbl-proimmes-models}

```{r spain replication proimmes models, echo = FALSE, warnings = FALSE}
#| label: tbl-proimmes-models

suppressWarnings(
  modelsummary(models_proimmes, output = "gt", 
             coef_map = c('treat1' = 'Treatment', 'imm_1' = 'Immigration', 
                          '(Intercept)' = 'Intercept'),
             gof_omit = "BIC|AIC|R2 Adj.|F|RMSE|Log.Lik.",  
             stars = c('*'=.1, "**"=.05, "***"=.01), 
             add_rows = tribble(~term,  ~Base,  ~Base, ~Base, ~Base, ~Base,
               "Weighted", "Yes", "No", "Yes", "No", "Yes", 
               "HC3 Robust SEs", "No", "No", "Yes", "Yes", "No",
               "Survey Robust SEs", "No", "No", "No", "No", "Yes"
             ))
)

```

Study 2 (Spain): Pro-Immigration Only Model (Table A9, Column 3) Sensitivity

:::

::: {#smtbl-noproimmes-models}

```{r spain replication noproimmes models, echo = FALSE, warnings = FALSE}
#| label: smtbl-noproimmes-models

suppressWarnings(
  modelsummary(models_noproimmes, output = "gt", 
             coef_map = c('treat1' = 'Treatment', 'imm_1' = 'Immigration', 
                          '(Intercept)' = 'Intercept'),
             gof_omit = "BIC|AIC|R2 Adj.|F|RMSE|Log.Lik.",  
             stars = c('*'=.1, "**"=.05, "***"=.01), 
             add_rows = tribble(~term,  ~Base,  ~Base, ~Base, ~Base, ~Base,
               "Weighted", "Yes", "No", "Yes", "No", "Yes", 
               "HC3 Robust SEs", "No", "No", "Yes", "Yes", "No",
               "Survey Robust SEs", "No", "No", "No", "No", "Yes"
            ))
)

```

Study 2 (Spain): Anti-Immigration Only Model (Table A9, Column 4) Sensitivity

:::

{{< pagebreak >}}

## Re-analyses of Western Values

In @smtbl-noweights-ancillary I re-analyse the results from Table A11 in the supplementary materials, which reports results for the ancillary and placebo outcomes. The only change here is to remove the weights, given that the original version of Table A11 uses robust standard errors. As column 2 shows, the results for Western liberal values disappears once weights are removed. 

::: {#smtbl-noweights-ancillary}

```{r spain replication no weights ancillary table, echo = FALSE, warnings = FALSE}
#| label: smtbl-noweights-ancillary

suppressWarnings(
modelsummary(mech, output = "gt", 
              coef_map = c('treat1' = 'Treatment', 'imm_1' = 'Immigration', 'treat1:imm_1' = 'Treatment x Immigration',
                          '(Intercept)' = 'Intercept'),
             gof_omit = "BIC|AIC|R2 Adj.|F|RMSE|Log.Lik.",  
             stars = c('*'=.1, "**"=.05, "***"=.01), 
             add_rows = tribble(~term,  ~Base,  ~Base, ~Base, ~Base, ~Base,  ~Base,
               "Weighted", "No", "No", "No", "No", "No", "No", 
               "HC3 Robust SEs", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes",
               "Survey Robust SEs", "No", "No", "No", "No", "No", "No"
             ))
)

```

Study 2 (Spain): Regression Results Without Weights for Ancillary and Placebo Outcomes (Re-Estimation of Table A11)

::: 

{{< pagebreak >}}

## Interaction effects by weight bin

In @smtbl-weight-hetfx-int I present the results of the interaction model estimated on three sub-samples -- those with low weights (under 0.01), those with high weights (over 3), and those mid weights (in-between). 

::: {#smtbl-weight-hetfx-int}

```{r weight bin hetfx_int table, echo = FALSE, warnings = FALSE}
#| label: smtbl-weight-hetfx-int

suppressWarnings(
modelsummary(models_hetfx_int, output = "gt",
            coef_map = c('treat1' = 'Treatment', 'imm_1' = 'Immigration', 'treat1:imm_1' = 'Treatment x Immigration', 
            '(Intercept)' = 'Intercept'),
            gof_omit = "BIC|AIC|R2 Adj.|F|RMSE|Log.Lik.",  
             stars = c('*'=.1, "**"=.05, "***"=.01), 
             add_rows = tribble(~term,  ~Base,  ~Base, ~Base,
               "Weighted", "No", "No", "No",  
               "HC3 Robust SEs", "Yes", "Yes", "Yes",
               "Survey Robust SEs", "No", "No", "No", 
             ))
)
```

Study 2 (Spain): Heterogeneous Interaction Effects By Weight Bin

:::

{{< pagebreak >}}

## Heterogeneous effects by age group

As noted in the main body of the comment, age is one of the primary factors that appears to differ between weight bins. A more detailed visualisation of the distribution of age across weight bins is presented in @smfig-age-by-bin. Essentially, every respondent who is over the age of 52 receives a high weight, those over the age of 35 but under 53 are much more likely to receive high weights than low weights, and those under 35 are much more likely to receive a low weight.  

::: {#smfig-age-by-bin}

```{r age plot by bin, echo = FALSE, warnings = FALSE}
#| label: smfig-age-by-bin

ggplot(spain, aes(x=age, y=weightbins)) +
  geom_jitter(height=.2, width=.2, alpha = .5, na.rm=TRUE) +
  ylab("Weight Bin") +
  xlab("Age") +
  theme_bw()
```

Study 2 (Spain): Age Distributions by Weight Bin

:::

One might imagine that age is correlated with immigration sentiment, such that older people are more predisposed to be anti-immigrant than younger people. This is not the case in the study 2 data. As shown in @smfig-imm1-by-age-es, there is no meaningful association between `age` and `imm_1`, and in fact the oldest respondents are often the most positively disposed toward immigration. 

```{r age and imm_1, echo = FALSE, include = FALSE, warnings = FALSE}
# standard error function:
stderr <- function(x, na.rm=FALSE) {
  if (na.rm) x <- na.omit(x)
  sqrt(var(x)/length(x))
}

# aggregate up imm_1 by age, Spain:
age_imm_es <- spain %>%
  filter(!is.na(age)) %>%
  group_by(age) %>%
  summarise(imm_1_mean = mean(imm_1, na.rm=TRUE),
            imm_1_se = stderr(imm_1, na.rm=TRUE), 
            count = n()) %>%
  ungroup()

```

::: {#smfig-imm1-by-age-es}

```{r age and imm_1 plot, echo = FALSE, include = TRUE, warnings = FALSE}
#| label: smfig-imm1-by-age-es

# plot:
ggplot(age_imm_es, aes(x = age, y = imm_1_mean)) +
  geom_point() + 
  geom_errorbar(aes(ymin = imm_1_mean - 1.96 * imm_1_se, ymax = imm_1_mean + 1.96 * imm_1_se), width = 0) +
  ylab("Average Immigration Attitude (Higher More Positive)")+
  xlab("Age") +
  theme_minimal() 

```

Study 2 (Spain): Age Is Unrelated to Immigration Sentiment

:::

A closer analysis is offered in @smtbl-hetfx-age, which shows the results of a simple linear regression of the outcome variable on the treatment indicator, subset into six age categories that were provided in the replication data. The treatment effect is only non-zero in the oldest age categories (45 years and older). The only statistically significant result is for the 45-54 category, and the point estimate is very large. While the point estimates for those in the two oldest categories are also quite large (albeit half the size of the 45-54 category), there are too few observations in this categories to say much with confidence. Notably it appears from this analysis that a small number of observations (around 150) likely drive the overall result presented in the paper. When these individuals are not up-weighted (and those who do not respond to the treatment are not down-weighted), the results unsurprisingly mostly attenuate. 

```{r hetfx by age analysis, echo = FALSE, warnings = FALSE}
# use the age categories already defined in the replication data:
models_hetfxage <- list(
  'Age < 25' = lm(support ~ treat, data=spain[spain$agecat==1,]),
  'Age 25-34' = lm(support ~ treat, data=spain[spain$agecat==2,]),
  'Age 35-44' = lm(support ~ treat, data=spain[spain$agecat==3,]),
  'Age 45-54' = lm(support ~ treat, data=spain[spain$agecat==4,]),
  'Age 55-64' = lm(support ~ treat, data=spain[spain$agecat==5,]),
  'Age >64' = lm(support ~ treat, data=spain[spain$agecat==6,])
)
```

::: {#smtbl-hetfx-age}

```{r hetfx by age table, echo = FALSE, warnings = FALSE}
#| label: smtbl-hetfx-age

modelsummary(models_hetfxage, output = "gt", robust = TRUE, 
             coef_map = c('treat1' = 'Treatment', '(Intercept)' = 'Intercept'), 
             gof_omit = "BIC|AIC|R2 Adj.|F|RMSE|Log.Lik.",  
             stars = c('*'=.1, "**"=.05, "***"=.01))
```

Study 2 (Spain): Heterogeneous Treatment Effects by Age Category

:::

It appears that the results in study 2 are largely on account of a combination of the weighting scheme targeting a small segment of the sample that appears to respond strongly to the stimulus. While this may be attributable to those individuals' characteristics (for example, that older people are theoretically more likely to respond to homonationalist appeals), it may well be down to pure chance. 

```{r imm_1 by age uk, echo = FALSE, warnings = FALSE}
# aggregate up imm_1 by age, UK:
age_imm_uk <- uk %>%
  filter(!is.na(age)) %>%
  group_by(age) %>%
  summarise(imm_1_mean = mean(imm_1, na.rm=TRUE),
            imm_1_se = stderr(imm_1, na.rm=TRUE), 
            count = n()) %>%
  ungroup()
```

::: {#smfig-imm1-by-age-uk}

```{r imm_1 by age uk, echo = FALSE, warnings = FALSE}
#| label: smfig-imm1-by-age-uk

# plot:
ggplot(age_imm_uk, aes(x = age, y = imm_1_mean)) +
  geom_point() + 
  geom_errorbar(aes(ymin = imm_1_mean - 1.96 * imm_1_se, ymax = imm_1_mean + 1.96 * imm_1_se), width = 0) +
  ylab("Average Immigration Attitude (Higher More Positive)")+
  xlab("Age") +
  theme_minimal()
```

Study 1 (UK): Age Is Unrelated to Immigration Sentiment

:::

Two things weigh against this deeper interpretation of these results. First, as was shown in @smfig-imm1-by-age-es, it is not the case that older voters are more likely to be more anti-immigration in the study 2 data. Second, for the UK (study 1) there is likewise no association between age and immigration sentiment, nor are there any coherent age heterogeneities. The relationship between age and immigration sentiment in the UK is shown in @smfig-imm1-by-age-uk, and the results of the binned regression analysis are shown in @smtbl-hetfx-age-uk. In this context, the youngest age category has a statistically significant positive treatment effect, while the second youngest has a statistically significant negative effect. The remaining categories are not statistically significant, and the point estimates bounce around between 0 and 0.1, with no clear pattern. 

```{r hetfx by age uk analysis, echo = FALSE, include = FALSE, warnings = FALSE}
# use the age categories already defined in the replication data:
models_hetfxage_uk <- list(
  'Age < 25' = lm(support ~ treat, data=uk[uk$agecat=="18-24",]),
  'Age 25-34' = lm(support ~ treat, data=uk[uk$agecat=="25-34",]),
  'Age 35-44' = lm(support ~ treat, data=uk[uk$agecat=="35-44",]),
  'Age 45-54' = lm(support ~ treat, data=uk[uk$agecat=="45-54",]),
  'Age 55-64' = lm(support ~ treat, data=uk[uk$agecat=="55-64",]),
  'Age >64' = lm(support ~ treat, data=uk[uk$agecat=="65+",])
)
```

::: {#smtbl-hetfx-age-uk}

```{r hetfx by age uk table, echo = FALSE, warnings = FALSE}
#| label: smtbl-hetfx-age-uk

modelsummary(models_hetfxage_uk, output = "gt", robust = TRUE, 
             coef_map = c('treat1' = 'Treatment', '(Intercept)' = 'Intercept'), 
             gof_omit = "BIC|AIC|R2 Adj.|F|RMSE|Log.Lik.",  
             stars = c('*'=.1, "**"=.05, "***"=.01))
```

Study 1 (UK): Heterogeneous Treatment Effects by Age Category

:::

{{< pagebreak >}}

## Corrected Figure 6 (No Weights)

In @smfig-6-correction-noweights I present the corrected (linear regression, robust standard errors) version of Figure 6 from the published paper, but with weights removed. 

```{r spain interaction no weights, echo = FALSE, include = FALSE, warnings = FALSE}
# run the linear model without weights
modelES_noweight <- lm(support ~ treat*imm_1, data=spain)
# use lm_robust for robust SEs that can be used by margins(). Set se_type = "HC3" to be consistent with summ(.,robust=TRUE). Point estimates are of course numerically identical bar rounding.
modelES_noweight_robust <- estimatr::lm_robust(support ~ treat*imm_1, data=spain, se_type = "HC3")

# generate code per replication materials, adding in the CI and making SE robust
predES <- interact_plot(modelES_noweight, pred = imm_1, modx = treat, interval = TRUE, robust = TRUE,
                        colors = colors)+
  labs(title="",
       y="Predicted support for\nLGBT+ education in schools",
       x="")+
  theme_minimal()+
  theme(legend.position = "none",
        axis.text.x =element_blank())+
  annotate(
    geom="text", x = 2.5, y = .65, size = 4, color = "#d11141", fontface=2,
    label = "Slope for\ntreated group")+
  annotate(
    geom = "curve", x =1.6, y = .6, xend = 1.2, yend = .44, 
    curvature = .4, arrow = arrow(length = unit(2, "mm")), colour="#d11141")+
  annotate(
    geom="text", x = 6, y = .4, size = 4, color = "#205C8A", fontface=2,
    label = "Slope for\ncontrol group")+
  annotate(
    geom = "curve", x =5, y = .4, xend = 2.6, yend =.31, 
    curvature = -.4, arrow = arrow(length = unit(2, "mm")), colour="#205C8A")

gg_df <-
  # move to robust object
  modelES_noweight_robust %>%
  margins(at = list(imm_1 = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treat1")

ameES<- ggplot(gg_df, aes(imm_1, AME)) +
  geom_point(colour="#d11141") +
  geom_line(colour="#d11141") +
  coord_cartesian(xlim = c(0, 10)) +
  # change to 95% ci:
  geom_errorbar(aes(ymax = (AME-SE*1.96), ymin = (AME+SE*1.96)), width = 0, colour="#d11141") +
  geom_hline(yintercept = 0, linetype = "dashed", colour="#205C8A") +
  xlab("Pre-treatment attitudes towards immigration")+
  ylab("Conditional ATE") +
  theme_minimal()

```

::: {#smfig-6-correction-noweights}

```{r spain interactions no weights plot, echo = FALSE, warnings = FALSE, fig.width=10, fig.height=6}
#| label: smfig-6-correction-noweights

predES/ameES+ 
  plot_annotation(title = 'Conditional average treatment effect: Study 2 (Spain)',
                  theme = theme(plot.title = element_text(size = 14, face="bold")))
```

Study 2 (Spain): Corrected Reproduction of Figure 6 Without Weighting

:::

{{< pagebreak >}}

# All Code Used In This Comment

This supplementary material was produced as a dynamic document using `Quarto`. All code run is included below.

```{r ref.label=knitr::all_labels()}
#| echo: true
#| eval: false
```
